classdef SO3FunMLS < SO3Fun
% a class representing a function on the rotation group
% Syntax
%   SO3F = SO3FunMLS(nodes,values)
%   SO3F = SO3FunMLS(nodes,values,N,__)
%   SO3F = SO3FunMLS(nodes,values,eps,__)
%   SO3F = SO3FunMLS(nodes,values,eps,w,__)
%
% Input
%  nodes  - @orientation,@rotation (interpolation points)
%  values - array of function values
%  N      - specified number of neighbors used for local interpolation
%  eps    - support radius of the weight function
%  w      - @function_handle (weight function)
%
% Output
%  SO3F - @SO3FunMLS
%
% Options
%  degree - the polynomial degree used for approximation
%
% Flags
%  centered     - only evaluate the basis near the pole if true
%  tangent      - use polynomials on the tangent space
%  hat          - use hat function as weight function
%  squared_hat  - use squared-hat function as weight function
%  indicator    - use indicator function as weight function
%

% TODO: transform into local interpolation-class where SO3FunMLS is a specific subclass

  properties
    nodes       = [];   % orientations where the function values are known
    values      = [];   % the corresponding values
    degree      = 3     % the polynomial degree used for approximation
    delta       = 0     % support radius of the weight function
    nn          = 0     % specified number of neighbors to use 
    w           = @(t)(max(1-t, 0).^4 .* (4*t+1)); % wendland weight function
    centered    = false % only evaluate the basis near the pole if true
    tangent     = false % use polynomials on the tangent space
    bandwidth   = getMTEXpref('maxSO3Bandwidth');
  end

  properties (Dependent)
    dim
    antipodal
    SLeft
    SRight
    isReal
  end

  methods
    
    function SO3F = SO3FunMLS(nodes, values, varargin)
    % initialize a SO(3)-valued function
    
      if nargin == 0, return; end
    
      % convert arbitrary SO3Fun to SO3FunHarmonic
      if isa(nodes,'function_handle') || isa(nodes,'SO3Fun')
        if nargin == 1, values=[]; end
        SO3F = SO3FunMLS.approximate(nodes,values,varargin{:});
        return
      end

      % TODO: uniqueData
      
      if isa(nodes,'rotation'), nodes = orientation(nodes); end
      SO3F.nodes = nodes; % preserve grid structure
      sz = [size(values),1];
      SO3F.values = reshape(values(:) , [length(nodes) , sz(find(cumprod(sz)==length(nodes), 1)+1:end)] );

      % set optional arguments
      SO3F.degree = get_option(varargin,'degree',3);
      
      % apply flags in the function arguments 
      SO3F.centered = get_option(varargin,'centered', false, 'logical');
      if get_option(varargin,'tangent', false, 'logical')
        SO3F.tangent = true;
        SO3F.centered = true;
      end

      % set the weight function 
      weightfun = get_option(varargin, 'weight', 'wendland', {'string','function_handle'});
      if (class(weightfun) == 'function_handle')
        SO3F.w = weightfun;
      else
        switch weightfun
          case 'hat';         SO3F.w = @(t)(max(1-t, 0));
          case 'squared hat'; SO3F.w = @(t)(max(1-t, 0).^2);
          case 'indicator';   SO3F.w = @(t)(t .* (t <= 1));
          case 'const';       SO3F.w = @(t)(t .* (t <= 1));
          case 'cos';         SO3F.w = @(t)((1+cos(pi*t))/2);
          case 'C1hat';       SO3F.w = @(t)((1-t.^2).^2);
          otherwise;          SO3F.w = @(t)(max(1-t, 0).^4 .* (4*t+1));
        end
      end

      % set delta or k if given
      if nargin > 2 && isnumeric(varargin{1})
        temp = varargin{1};
      else
        temp = 2 * SO3F.dim;
      end
      % if the input is a whole number, assume that nn is specified
      if (floor(temp) == temp)
        if (temp < SO3F.dim)
          warning('The specified number of neighbors nn was less than the dimension dim. nn has been set to 2 * dim.');
          SO3F.nn = 2*SO3F.dim;
        else 
          SO3F.nn = temp;
        end
        SO3F.delta = guess_delta(SO3F);
      else
        SO3F.nn = 2 * SO3F.dim;
        SO3F.delta = temp;
      end
      
    end

    function dimension = get.dim(SO3F)
      if (SO3F.degree == 0)
        dimension = 1;
        return;
      end
      dimension = nchoosek(SO3F.degree + 3, 3);
    end

    function d = guess_delta(SO3F)
      % for N nodes on one hemisphere, the expected number of nodes in a
      % spherical cap of angular radius phi is
      %         N * 2/pi * (phi - sin(phi) * cos(phi))
      % choose delta such that the expected number of neighbors is 2*sF.dim
      syms phi;
      d = double(vpasolve(phi-sin(phi)*cos(phi) - pi*SO3F.dim/numel(SO3F.nodes)));
      % the quaterion distance is twice the spherical distance
      d = 2 * d;
    end

    % if only delta is specified, guess nn for this delta
    function nn = guess_nn(SO3F, varargin)
      testnodes = equispacedSO3Grid(SO3F.nodes.CS, 'points', 1000);
      ind = SO3F.nodes.find(testnodes, SO3F.delta);
      if (numel(varargin) == 0)
        nn = ceil(mean(sum(ind, 2)));
        return;
      end

      if (varargin{1} == "min")
        nn = floor(min(sum(ind, 2)));
      elseif (varargin{1} == "max")
        nn = ceil(max(sum(ind, 2)));
      end
    end

    function SO3F = set.SRight(SO3F,S)
      SO3F.nodes.CS = S;
    end

    function S = get.SRight(SO3F)
      try
        S = SO3F.nodes.CS;
      catch
        S = specimenSymmetry.default;
      end
    end

    function SO3F = set.SLeft(SO3F,S)
      SO3F.nodes.SS = S;
    end

    function S = get.SLeft(SO3F)
      try
        S = SO3F.nodes.SS;
      catch
        S = specimenSymmetry.default;
      end
    end

    function SO3F = set.antipodal(SO3F,antipodal)
      SO3F.nodes.antipodal = antipodal;
    end

    function antipodal = get.antipodal(SO3F)
      try
        antipodal = SO3F.nodes.antipodal;
      catch
        antipodal = false;
      end
    end

    function out = get.isReal(f)
      out = isreal(f.values);
    end
  
    function F = set.isReal(F,value)
      if ~value, return; end
      F.values = real(F.values);
    end


    function n = numArgumentsFromSubscript(varargin)
      n = 0;
    end

  end

  methods (Static = true)
    SO3F = interpolate(varargin);
    SO3F = approximate(f, varargin);
    SO3F = example(varargin)
  end
  
end
